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Abstract 

One of the objectives of a Delaunay mesh refinement algorithm is to produce meshes with 
tetrahedral elements having a bounded aspect ratio, which is the ratio between the radius 
of the circumscribing and inscribing spheres. The refinement is carried out by inserting 
additional Steiner vertices inside the circumsphere of a poor-quality tetrahedron (to remove 
the tetrahedron) at a sufficient distance from existing vertices to guarantee the termination 
and size optimality of the algorithm. This technique eliminates tetrahedra whose ratio 
of the radius of the circumscribing sphere and the shortest side, the radius-edge ratio, 
is large. Slivers, almost-planar tetrahedra, which have a small radius-edge ratio, but a 
large aspect ratio, are avoided by small random perturbations of the Steiner vertices to 
improve the aspect ratio. Additionally, geometric constraints, called “petals”, have been 
shown to produce smaller high-quality meshes in 2D Delaunay refinement algorithms. In 
this paper, we develop a deterministic nondifferentiable optimization routine to place the 
Steiner vertex inside geometrical constraints that we call “snow globes” for 3D Delaunay 
refinement. We explore why the geometrical constraints and an ordering on processing of 
poor-quality tetrahedra result in smaller meshes. The stricter analysis provides an improved 
constant associated with the size optimality of the generated meshes. Aided by the analysis, 
we present a modified algorithm to handle boundary encroachment. The final algorithm 
behaves like an advancing-front algorithms that are commonly used for quadrilateral and 
hexahedral meshing. 


‘The work of the author was supported in part by the NIH/NIGMS Center for Integrative Biomedical Com¬ 
puting grant 2P41 RR0112553-12 and a grant from ExxonMobil. The author would also like to thank Christine 
Pickett, an editor at the University of Utah, for proofreading and finding numerous typos in the paper. 
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1 Introduction 


Delaunay refinement is usually carried out by eliminating one poor-quality tetrahedron at a 
time. Jumpani and Ungor [B] extended the work to eliminate tetrahedra associated with a 
single short edge in the mesh by simultaneous insertion of multiple vertices around the edge. In 
this paper, we extend the algorithm to simultaneously eliminate tetrahedra all over the mesh by 
multiple Steiner vertex insertions. In addition, we replace the exisiting randomized algorithm 
to avoid slivers, which are almost-planar tetrahedra having small dihedral angles and a large 
aspect ratio, with a combinatorial algorithm for the optimization of nondifferentiable objective 
functions. 

Delaunay refinement is carried out by placing additional Steiner vertices inside the circum- 
spheres of poor-quality elements to obtain a higher-quality mesh EEiEiiiniiiaiiaE]. The 
radius-edge ratio, the ratio between the circumscribing circle and the shortest edge of a tetra¬ 
hedra, may be used to measure the quality of an element. If the ratio is large, the element 
is “split” by adding its circumcenter to the mesh (and re-tetrahedralizing the domain). If the 
ratio is greater than 1, the new vertex introduces edges that are at least as large as the ex¬ 
isting shortest edge in the tetrahedron. Thus, the termination of the algorithm can be proved 
by the packing argument. To account for the elements near the boundary (where constraints 
may prevent the insertion of new vertices), the ratio has to be greater than 2-\/2/3 [15]. If a 
Steiner vertex results in a large sliver, the vertex is added to the mesh because a subsequent 
insertion of a vertex at the circumcenter of the sliver eliminates the sliver and introduces a long 
enough edge that still guarantees the termination of the algorithm. For small slivers, however, 
the new vertex is perturbed enough such that it introduces a long-enough edge and its location 
results in a better-quality tetrahedron. An analysis of the quality of tetrahedra produced in 
such adversarial cases provides the theoretical bounds on the mesh quality IZIEIIH]. 

The perturbation of the vertex location may be formulated as a nondifferentiable optimiza¬ 
tion problem. In our recent work m , we developed a combinatorial algorithm for such problems. 
We describe this algorithm and some background on mesh generation in Section 2. Some 2D 
off-center vertex insertion algorithms for Delaunay refinement place vertices to the extent pos¬ 
sible from other vertices, but within certain geometric constraints (called “petals”) [1711181 H]. 
These algorithms place only one vertex at a time. We describe these algorithms in Section 3. 
In our algorithm in Section 4, we avoid constructing small slivers as far as possible by placing 
vertices using our combinatorial optimization algorithm. We also provide an iterative frame¬ 
work so that multiple vertices can be added simultaneously to further improve the mesh quality. 
The 2D geometric constraints (petals) get translated to “snow globes” in our 3D algorithm. In 
Section 5, we carry out geometric analyses to show that our algorithm terminates and produces 
size-optimal meshes. For single-vertex insertion, the proofs were provided in earlier work, but 
we carry out a stricter analysis in which we account for the fact that we process poor-quality 
tetrahedron with shorter edges first. Thus, we obtain an improved constant of proportionality 
associated with the size optimality of the mesh. This analysis is an extension of Si’s work m- 
This analysis also helped us develop a modified algorithm to handle boundary facets and edges. 
For multiple-vertex insertion, we prove that small perturbations in vertex positions due to our 
combinatorial optimization algorithm still satisfy the conditions necessary to prove the termi¬ 
nation and size optimality. Possible future research directions are discussed in Section 6. 

2 Background 

We present some notations and background on mesh quality metrics, Delaunay refinement, and 
optimization of nondifferentiable functions. 
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(a) “top” view 


(b) “front” view 


(c) potential slivers 


Figure 1: The first two figures provide the potential locations of the fourth vertex such that 
it forms a sliver with three other vertices. The cross-section of the “forbidden” region has 
an hourglass-like shape. The third figure depicts that there may be only a finite number of 
potential small slivers for an almost-good mesh inside the circumsphere of a tetrahedron. 

2.1 Mesh Quality 

We begin by defining the following quality metrics and other terms commonly used in mesh 
generation algorithms. 

Definition 2.1. Radius-Edge Ratio: The radius-edge ratio, p, of a simplicial element 
is the ratio of the radius of its circumsphere and its shortest edge. 

Definition 2.2. Volume-Edge Ratio: The volume-edge ratio, a, of a simplicial ele¬ 
ment is the ratio of its volume and the cube of its shortest edge. 

Definition 2.3. Almost-Good Mesh: A mesh is /9*-almost-good if the radius-edge 
ratio of all its elements is less than a constant p*. 

Definition 2.4. Sliver: A sliver is a tetrahedron whose radius-edge ratio is less than 
a constant p*, and the volume-edge ratio is less than a constant a*. 

Definition 2.5. Local Feature Size: The local feature size at a point x is the smallest 
disk/sphere that can be placed at x such that at least two nonincident features are 
enclosed by the disk/sphere. 

Given a triangle qrs, the region of potential locations of an off-plane vertex p such that it 
forms a sliver pqrs is inside the region shown in Fig. [H called the “forbidden” region. If I is the 
shortest side of Aqrs, the forbidden region is bounded by spheres whose radius is p*l and two 
planes parallel to the plane of Aqrs at approriate distances such that the volume-edge ratio is 
below a*. Li [7] showed that there are at most a finite number of smal0 slivers that may be 
formed when a vertex is being inserted near the circumcenter to refine an almost-good mesh. 

The size optimality of an algorithm is typically shown using the local feature size, lfs{x) at 
any point x in the domain. The size-optimality of a mesh is proved by showing that the length 
of each edge has a lower bound proportional to the local feature size. This lower bound shows 
that the mesh obtained by an algorithm is, at most, a constant times the size of the smallest 
possible mesh of the desired quality. 

2.2 Nondifferentiable Optimization 

We provide a brief description of the first-order necessary conditions for the existence of a local 
optimum of nondifferentiable function at a point. Consider a composite function f{x) that is 
defined as the minimum over a set of n differentiable functions fi{x) at any point x, Vi G [1, n]. 
Clearly, / is a nondifferentiable function, and we wish to find x* that maximizes f{x). A vector 
g is dehned as a subgradient of / at x if 3 an e neighborhood such that F{x -be) — /(x) > ^£. 

^If the circumradius of a sliver is below a certain threshold, it is called a small sliver. 
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If only one function fi, i G [1, n] defines the value of f{x) at x, the gradient of /* at x is also 
the subgradient of f{x) at x. If multiple function^ define the value of f{x) at a point, any 
convex combination of the gradients of the active set functions is also a subgradient of / at x. 
The set of all subgradients is called the subdifferential, df, of / at x. 

For unconstrained optimization, a necessary condition for a local minimum is 0 E df, 
i.e., some convex combinations of the gradients of the active set functions should vanish. For 
constrained optimization, the first-order necessary condition includes the gradients of the active 
constraints in the convex combination. Formally, let hi(x) < 0 be a set of k constraints, 
i E [1, k], that the location x* has to satisfy. At a local minimum, the Kahn-Karush-Tucker 
(KKT) conditions have to be satisfied, i.e., hi{x*) < 0, A* > 0, X*hi{x*) = 0, and 0 E df{x*) — 

For a d-dimensional function, if m > d-|-I subgradients have a vanishing convex combination, 
there exists a set of d -I- 1 subgradients (among the m subgradients) that also have a vanishing 
convex combination. Thus, it is sufficient to consider all possible combinations of d -|- 1 or 
fewer subgradients to verify if a convex combination vanished. In our algorithm, we consider 
all possible combinations of d -|- 1 or fewer constituent functions, fi{x), and constraints, hi{x). 
For each combination, we find locations where the constituent functions are equal and intersect 
all the constraints. The locations are, typically, determined by solving a system of d -|- 1 
polynomial equations (or nonlinear equations) with d -|- 1 variables (d dimensions and the value 
of the constituent function), which has a finite number of solutions. We consider all solutions 
for each combination and return the optimal solution. See m for an example involving the 
maximization of the minimum angle. 

3 Related Work 

In this section, we highlight prior research in Delaunay refinement. 

3.1 Delaunay Refinement 

In order to ensure that the algorithm terminates, the vertex has to be inserted at least a finite 
distance away from all other existing vertices. If the radius-edge ratio is p, the circumcenter is 
at least a distance pi from all other vertices of the mesh, where I is the length of the shortest side 
of the triangle. In fact, any point inside the circle/sphere with the circumcenter as the center 
and having a radius of {p — 1)1 (for p > 1) is, at least, a distance I from other vertices. This 
sphere/circle is called a “picking” region. In our algorithm, we choose a parameter 1 < a < p* 
such that the radius of the picking region is {p — a)l. This ensure that every new edge is longer 
than the current shortest edge in the mesh. It was shown in [5] that inserting any point in the 
picking region is sufficient to obtain a size-optimal mesh (with the minimum angle greater than 
arcsinl/\/5 for 2D refinement). In 3D, a larger p and a smaller a results in a bigger sphere 
where there is more room for the Steiner vertex to avoid creating slivers, but a larger p also 
implies a poorer-quality tetrahedron and a smaller a may a bigger mesh with shorter edges. 
Thus, the radius-edge ratio and the volume-edge ratio is chosen so that the aspect ratio (or the 
dihedral angle) of the resulting tetrahedra should be maximized. Li [7] analyzes this trade-off. 

The algorithms above do not specify the exact location where the Steiner vertex should be 
placed so that the resulting mesh has fewer vertices and elements than the naive circumcenter 
insertion. In 2D, the answer to that question was provided by Ungor and Erten [numiu. In 
the first paper Ungor considers the user-defined minimum angle specification and carries 
out a targeted placement of vertices, which results in a triangulation whose minimum angle is 

^They are also called the active set functions. 

the origin is inside the convex hull of the vectors that represent the subgradients, a convex combination 
vanishes. 
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(a) first vertex insertion (b) second vertex insertion (c) simultaneous insertion 


Figure 2: The insertions of two vertices when they are constrained to be inside the “petals” 
constructed on the shortest sides of two poor-quality triangles. On the left (two figures), one 
vertex is inserted at a time. The first vertex is inserted on the circumference, and the second 
vertex is inserted at the circumcenter of the new triangle. On the right, both vertices are 
inserted simultaneously by maximizing the minimum distance. Notice that the triangulation 
on the right has lexicographically better angles. 


the user-specified value. The vertex is placed on the perpendicular bisector of the shortest edge 
of a poor-quality triangle such that the angle subtended by the edge on the point is equal to the 
required minimum angle. If the point of intersection is too close to other vertices, the algorithm 
reverts to the circumcenter insertion technique. This algorithm produces meshes with nearly 
50% fewer vertices. 

In Erten and Ungor’s paper [1], they place the Steiner vertex within a “petal” (see Fig. ED, 
which is the locus of the points at which the shortest edge subtends the required minimum angle, 
at a location that maximizes the minimum distance from other vertices. Since this choice of the 
vertex location is within the constraints of the previous algorithm [5], it can be shown that the 
technique terminates and produces size optimal meshes with a minimum angle of arcsin 1 / ^/5. 
In Section 5, we show why placing vertices inside the petals associated with the shortest edges 
first results in meshes having few vertices and elements. This algorithm produces meshes with 
nearly 50% fewer vertices than the previous algorithm. 

The algorithm in the first paper m was extended to 3D [B] by placing the vertices on the 
perpendicularly bisecting plane of the shortest side. Our algorithm extends the idea of the 
second algorithm to 3D by defining a “snow globe”, which is analogous to a petal, in which we 
place Steiner vertices as far away as possible from existing vertices while avoiding the forbidden 
regions defined in Section 2. 

3.2 Simultaneous Multiple-Vertex Insertion 

In our recent research [12], we extended Erten and Ungor’s algorithm [3| such that multiple 
vertices may be simultaneously placed by the algorithm. Their insertion algorithm maximizes 
the minimum distance to existing vertices in a greedy manner, i.e., each vertex is placed at 
such a location by considering existing vertices in the mesh, but not the vertices that may be 
placed later. In our extension, we consider simultaneously placing vertices inside some of the 
current petals. We begin with the greedy placement strategy as in the previous algorithm and 
iteratively maximize the minimum distance. Fig. [2] provides the motivation for this algorithm 
by providing a comparison of the greedy strategy with our strategy to place the vertices. It can 
be seen that our strategy provides lexicographically better angles than the greedy strategy. 

Lemma 3.1. Given a 2D input vertex set, one of the Voronoi vertices or one of the points of 
intersection of the Voronoi edges and the convex hull maximizes the minimum distance to any 
vertex of the input set if the loeation is eonstrained to he within the convex hull of the input 
vertex set. 

Proof. In order to find the location that maximizes the minimum distance, consider the gradients 
of the distance function from each of the input vertices. The function is nondifferentialble at 
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(a) local maximum (b) saddle point (c) constrained maximum (d) “stationary” points 


Figure 3: Stationary points to determine the minimum distance between any two points for 
inserting one additional vertex. The blue, green, and red vertices are the input vertices, the 
black vertices are optimal locations, and the arrows indicate the gradients. The gradients are 
directed such that their convex combination may vanish. In the last two figures, the vertex is 
constrained to be in the gray regions. In the last figure, although the gradients are suitably 
directed at the two marked points, only one of the points (at the top) is a local maximum. 


the Voronoi edges (and vertices) as distance functions from multiple vertices interact at these 
locations. As described in Section 2, a local maximum occurs when one, two, or three gradients 
have a vanishing convex combination. Since the distance function for a single vertex does not 
have a zero gradient, we consider the other two cases (see Fig. [3|). When two gradients have 
a vanishing convex combination, they are antiparallel, which happens at the midpoint of two 
vertices. This point is a saddle point if it is inside the convex hull, but it is a local maxima if it 
is on the edge of the convex hull. It is easy to show that three gradients have a vanishing convex 
combination at all the Voronoi vertices. Since we have considered all possible local mamima, 
one of these locations maximizes the minimum distance. □ 

In Erten and Ungor’s algorithm, additional potential locations of a local maximum are 
considered on the boundary of a petal as KKT conditions may be satisfied at these locations. 
In our extension |12] , we iteratively maximize the minimum distance for each new Steiner vertex 
(and constrain them to be within their respective petals) by allowing only one vertex to move 
in each step and maximizing its distance from all other vertices. We retriangulate the domain 
when we are satisfied with the vertex distribution. The iterative algorithm is similar to Lloyd’s 
algorithm to obtain centroidal Voronoi tessellation. The algorithm converges to a stationary 
point because we always maximize the minimum distance in each step. This algorithm can be 
extended to 3D in a straightforward manner. 


4 Deterministic Algorithms for Sliver-Pree Delaunay Refine¬ 
ment 

4.1 Preliminaries 

The Snow Globe For 2D refinement, petals were defined on the 

shortest side of a poor-quality triangle. For 3D refinement, we define 

the “smallest” facet as one of the two triangular faces, t, that contains 

the shortest side. It, of a poor-quality tetrahedron T and the next 

shortest side that is adjacent to If the circumradius of t is T), 

its radius-edge ratio is pt = Ytjlt- If Pt < P*, the radius-edge ratio 

threshold for an almost-good mesh, a snow globe is defined by the 

sphere of radius p*lt circumscribing the facet t as in Fig. 01 Note that 

there are two such spheres, and we consider the sphere with the larger 

part on the same side of t as the fourth vertex of the tetrahedron. If 
_ Figure 4: A snow globe 

^This reason for this choice becomes clear is Section 5. 
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the radius-edge ratio of all facets of a tetrahedron are too large, a 
snow globe may not be defined. In Section 5, we will elaborate of 
what needs to be done under such circumstances. 

Additional Geometric Constraints We know that the inserted Steiner point has to be at 
a distance alt from the vertices of facet t. Thus, a vertex is constrained to be within the picking 
region described in Section 3. In addition, the forbidden regions shown in Fig. [1] are also added 
as constraints to avoid each potential small sliver. 

Objective Function Formulation For the distance maximization-based routine below, the 
objective function is simply the maximization of the minimum Eucleadian distance between the 
Steiner vertices and the existing vertices. We solve systems of four four-variable polynomials, 
where the minimum distance is itself another variable to be computed. 

4.2 Single-Vertex Delaunay Refinement Algorithm 

The distance maximization-based routine is quite simple. As in |Tj, we place our vertices as 
far as possible from existing vertices within the snow globe and the picking region, but outside 
the forbidden region. Our combinatorial optimization algorithm may be used to maximize the 
distance under these constraints. If an inserted vertex encroached npon a boundary facet or an 
edge, all vertices with the diametral circle/sphere are deleted and the midpoint/circumcenter 
of the edge/facet is added to the mesh. The readers are directed to prior research [TTl 141 [5l fT5] 
to learn more about handling such cases. It is important to add vertices in the snow globes 
associated with shortest edges first to get small meshes. In the next section, we provide a 
technique to enforce this even for boundary-encroaching Steiner vertices and boundary elements. 
A disadvantage of using the distance-based technique is that the definition of the sliver is 
essential to avoid creating the slivers; the definition is explicitly used as a constraint. The angle 
maximization-based routine provided in the appendix does not have this disadvantage, but it 
may produce larger, better-quality meshes. 

4.3 Multiple-Vertex Delaunay Refinement Algorithm 

For multiple vertices, we use an iterative technique for the distance maximization routine. We 
consider several snow globes and maximize the minimum distance one by one iteratively. After 
each such vertex movement, the forbidden regions are updated. We want to ensure that poor- 
quality tetrahedra with shortest sides are processed first. Therefore, we follow an incremental 
vertex placement strategy, where we place an addition vertex and redistribute them until we 
run out of the current set of snow globes or the minimum distance is less than the shortest edge 
of the tetrahedra responsible for the current or next snow globe. The algorithm is presented 
Algorithm 1. Note that the a perimeter-based distance mentioned in the algorithm is the 
not just the picking region constraint, but also includes the distance from “candidate” Steiner 
vertices that are to be added in this step. 

5 Proofs and Guarantees 

The analysis of the algorithm in this section provides the reasoning for our choice of the ordering 
of poor-quality tetrahedra (and the corresponding snow globes) based on the shortest sides. The 

vertex encroaches upon a boundary facet or a boundary edge if it is inside the diametral “lemon” |14| of 
the facet or diametral “lens” [5] of the edge. A different definition is provided in in which a vertex encroaches 
upon a boundary edge/facet if it is invisible from the interior of the triangle/tetrahedron due to the boundary 
edge/facet. Either definition may be used. 
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Data: A Delaunay tetrahedralization of the domain with poor-quality tetrahedra. 

Result: A new triangulation of the domain with some of the poor-quality tetrahedra 
replaced with better-quality tetrahedra. 

compute the list of all poor-quality tetrahedra and the associated snow globes; 
while any of the snow globes are empty and the minimum distance among the next set of 
Steiner vertices is greater than the shortest edge associated with the next snow globe do 
insert a new Steiner vertex in the snow globe associated with the next shortest edge; 
move the Steiner vertices within their petals and maintain the minimum a 
parameter-based distance to locations that maximize the minimum distance outside 
the forbidden region; 

if any vertex encroaches upon an input facet/edge then 

delete all free vertices inside the diametral sphere/circle of all such input edges; 
insert the circumcenter/midpoint of the facet/edge; 

maximize the minimum distance among the remaining Steiner vertices within 
their snow globes and picking regions, but not encroaching upon any boundary 
facet/edge; 

break from the while loop; 

end 

end 

retetrahedralize the domain; 

Algorithm 1: An algorithm to simultaneously insert Steiner vertices for Delaunay refinement 


analysis also paves the way for handling vertex encroachment on the boundary facets and edges 
in a slightly different manner. The main difference between our analysis and that of Si [16] 
is that we consider the ordering and also constrain the new vertex to be within a certain 
distance dictated by the dimensions of the snow globes. This feature makes this technique an 
advancing front technique. Important: This analysis assumes that the length of the shortest 
edge adjacent to vertex v, denoted by r^, is Isf(u) or less for all vertices in the initial piecewise 
linear complex (PLC). We discuss the extension to a general case in the appendix. We also 
provide some additional analysis and discussion about our algorithm in the appendix. 

5.1 Termination 

For the single-vertex insertion algorithms, the proof of termination is identical to Foteinos et 
al.’s [5] work because we place the vertices within the geometric constraints defined by their 
algorithm. 

Lemma 5.1. The vertex insertion algorithm terminates for appropriate p* and a* values that 
are prescribed in the analysis by Li n- 

Proof. For p* > 2, the maximization of the minimum distance inside the snow globe yields 
an edge length of at least the shortest side of the poor-quality tetrahedra. Since only a finite 
number of small slivers are possible, there exists some a* for a point inside the picking region 
and the snow globe, but outside the forbidden region (corresponding to the a* and p*). Note 
that since we consider not only the picking region but also the intersection of the picking region 
with the snow globe, our bounds may be slightly different. We do not, however, analyze the 
bounds because it does not provide any additional insight into the problem, and the theoretical 
bounds are too small when compared with prior practical results. □ 

Lemma 5.2. The multiple-vertex insertion algorithm terminates. 




Proof. Every vertex movement in the distance-based algorithm increases the minimum distance 
between the vertex of interest and some other vertex. Thus, the minimum distance can only 
increase at every step. □ 

5.2 Size Optimality 

We will prove the size optimality for the 2D single-vertex insertion algorithm. Based on the 
analysis, we make suitable modifications to construct a 3D size-optimal algorithm. Then, the 
size optimality of the 3D algorithm is self-evident. We assume that the desired radius-edge 
ratio is p*, and the desired volume-edge ratio is a*. The idea behind the analysis is to show 
that when we order the insertion for petals associated with smaller edges first, the vertices are 
inserted in an advancing front manner. Thus, new vertices are inserted relatively close to the 
existing set of good-quality triangles. The new poor-quality triangles have longer edges, and 
they will insert vertices slightly further away. In this process, the lengths of the new shortest 
edges of poor-quality elements progressively increase. This progressive increase is absent when 
poor-quality triangles are ordered in a different manner or when the vertices are inserted away 
from the petals. We will show this result assuming that boundary encroachments do not play a 
role. Then, we make suitable modifications to our algorithm when boundary edges/facets come 
into play. 

Lemma 5.3. There are constants Ci and C 2 such that CiZmin < ^new < C 2 lmin, where Zmin is 
the length of the shortest side of the poor-quality triangle in consideration and Zmin is the length 
of the shortest new edge adjacent to the current shortest side after the insertion of the vertex 
in the petal. 

Proof. Since we enforce the picking region constraint, Ci = a. The size of the petal is a 
function of the shortest side of the mesh. Thus, the Steiner vertex may be inserted within a 
certain distance from the vertices on the shortest edge. We will call C 2 as /3 in the rest of the 
paper. If a snow globe exists, this lemma extends to 3D as well. □ 

As mentioned in the beginning of this section, we assume that /min is also the minimum 
local sizing field in the domain. For the sake of convenience, it is easy to view the vertex 
insertion procedure as taking place in stages. In the first stage, we insert vertices with the 
petals associated with poor-quality triangles of the shortest sides with lengths from /min to 
d/ min , where /min is the shortest edge in the mesh. In the second stage, the lengths of the 
shortest edges are from a/min to fPlmm (due to the lemma above). In the stage, the shortest 
sides have lengths from (^"'“^/min to /3^/min- 

Lemma 5.4. At each stage of the algorithm, the locations of the inserted Steiner vertices have 
constraints on the local sizing field at their respective locations. 

Proof. The local sizing field, Ifs(x) is a Lipschitz function, i.e., |lfs(x) — lfs(y)| < ||x — yHij- 
Thus, in the first state, if a Steiner vertex is inserted at point p, Isf(p) < /min + /3/min- In the 
n stage, Isf(p) < /min T fdlram T /? /min T /d^/min T ••• + /3 /min — ^ j^_i /min- D 

We define a protected region around each vertex as follows. It is easy to see that no new 
Steiner vertex may be added inside the protected region. 

Definition 5.1. Protected Region: A protected region is a circle/sphere centered at 
every vertex in the mesh having the radius of the shortest side of a poor-quality 
tetrahedron among all such tetrahedra at any stage of refinement. 

Lemma 5.5. The protected region grows exponentially after each stage. 
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Proof. The length of the shortest edge of a poor-quality tetrahedron in the first stage is /min- 
In the second stage, the length becomes a/min- In the the length becomes Thus, 

the protected region grows exponentially. □ 

Theorem 5.6. The single-vertex insertion algorithm produces size-optimal meshes. 

Proof. In order to prove this statement, we will use Si ’s m technique. Consider a vertex uq 
and the set of petals into which it was insertec^. Let us define the parent of uq as the closest 
vertex vi among the vertices to which the petals belong. Let the edge associated with the 
petal be of length Iq. We know that |uoUi| = XqIq, where a < Xq < (3. Let the parent of 
vi be V 2 and the parent of V 2 be U 3 and so on. Let the corresponding lengths of the shortest 
sides be h, I 2 , and so on. Let the corresponding ratio of the lengths \vkVk+i\/lk be A^. If such 
a sequence ends at Vm, where Vm is the input vertex, the local sizing field at vq is bounded 
from above by lm{l + Am + AmAm-i + AmAm-iAm -2 + AmAm-i.-.Ao). The length Iq is equal to 
lmXmXm-iXm-2---X0 by definition. The ratio of the local sizing field to the length Iq is given by 

1 1 1 1 ^ 1111^1 
AmAm-l.-.Ao Am-lAm- 2 -.-Ao Am- 2 --.Ao Aq ~ q, — 1 — (I/q.) 

Since the ratio is bounded and the protected region grows, the vertex closest to vq can only be 
placed at a distance of at least oIq. Thus, the ratio of the local sizing field at a vertex to the 
shortest edge edjecent to the vertex is given by 

1 1 _ 1 
1 — ( 1 /a) a a — 1 


, and we obtain a size-optimal mesh. □ 

The constant of proportionality associated with our algorithm is better than the one derived 
by Si dS] as we have considered the ordering of the poor-quality elements. Note that the analysis 
was carried out for 2D refinement without considering boundary encroachment. Below, we 
explain the necessary modifications for the 3D refinement algorithm and for handling boundary 
encroachment. 

5.2.1 3D Refinement 

For the 3D Delaunay refinement algorithm, when the snow globe is present for a triangular 
facet, the edge length-based conditions above hold. When the snow globe is not present, we 
have to enforce the conditions. Consider a poor-quality tetrahedron that does not have a snow 
globe. Consider its facet with the shortest edge and the next shortest edge connected to it. 
Construct a petal associated with the triangle for some radius-edge ratio p < p*, and then 
construct the spindle torus costructed by rotating the petal about the shortest edge. Add a 
vertex inside the spindle torus and the a parameter-based picking region, retetrahedralize the 
domain, and construct the snow globe on the new facet if necessary. 

Lemma 5.7. The intersection of the spindle torus and the picking region is non empty. 

Proof. Consider the plane formed by the shortest edge of the tetrahedron and its circumcenter. 
The projection of the spindle torus on the plane is a petal, and the projections of the circumcircle 
and picking region on the plane are two concentric circle. For p > \/2 and 1 < a < p, it is easy 
to show that the intersection of the petal and picking circle is nonempty. Thus, the intersection 
of the spindle torus and the picking region is nonempty. □ 

®A vertex may be part of multiple petals by chance. 
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V me 


(a) boundary edge 



(b) boundary facet 


Figure 5: If boundary elements are encroached upon, earlier methods placed (red) vertices on the 
midpoint of an edge or the circumcenter of a facet and deleted all vertices inside the diametral 
circle/sphere of the edges/facets. In our modification, we place the vertices (blue) close to the 
“advancing front” (but sufficiently far away from other vertices to guarantee termination and 
size optimality of the final mesh) and delete vertices inside the smaller circles/spheres as shown. 
On the left, the poor-quality 2D element is shown. On the right, only the boundary facet is 
shown. 


5.2.2 Handling Bonndary Encroachment 

Consider the queue that dictates the order of the petals/snow globe into which a Steiner vertex is 
inserted based on the shortest side associated with each petal/snow globe. If a potential Steiner 
vertex in a petal/snow globe encroaches upon a boundary edge or facet, the petal should be 
ordered in the queue as if its shortest side is either (a) its own shortest side Zmin or (b) half the 
length of the boundary edge or the radius of the circumcircle of boundary facet /mid divided by 
a, i.e., Iraid/cn-, whichever is shorter. Let us call this the effective length, /eS- The idea is to get 
rid of poor-quality triangles with the shortest edges first. Since such petals or snow globes can 
bring about poor-quality triangles with shorter edges, the petals are appropriately ordered. 

If a Steiner vertex does encroach upon the boundary edge/facet, and if there exists a point 
inside the petal/snow globe and the picking region, place the vertex at such a point. If not, first, 
we choose one of the vertices of the encroached boundary edge/facet that has the lowest local 
sizing field. From that vertex v, we choose a point m on the encroached boundary edge/facet 
such that \vm\ = 'jles, where a < 7 < /3 and \vm\ < /mid (see Fig.[S]). If we are handling a facet, 
we choose m on the line joining v and the circumcenter of the facet. Second, add m to the mesh, 
and delete all vertices within the circle/sphere formed with m as the center and \vm\ as the 
radius. Finally, retetrahedralize the domain, which creates new triangles/tetrahedra having the 
shortest side of length, at least, y/efr- Moreover, the vertex m is added close to the “advancing 
front”. Thus, our size optimality results still hold when boundary edges/facets are considered 
in the above lemmas and theorem. Notice that when /gg = /mid) the algorithm is equivalent to 
the prior Delaunay refinement algorithms. 

5.2.3 Size Optimality for the Multiple-Vertex Insertion Algorithm 

Theorem 5.8. The multiple-vertex insertion algorithm terminates with size-optimal meshes. 

Proof. In Algorithm 1, we simultaneously add multiple vertices only if the minimum distance 
is greater than a times the shortest edge associated with the next snow globe. Thus, the 
constraints on the minimum and the maximum distance from the current “front” are still 
maintained. We enforce the ordering on the snow globe explicitly through the suitable condition 
in the while loop in Algorithm 1 . Therefore, all of the lemmas above hold, and a size-optimal 
mesh is generated. □ 
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6 Future Work 


We have extended a 2D Delaunay refinement algorithm to 3D (petals to snow globes; deter¬ 
ministic maximization; multiple vertex insertion) and consolidated 3D algorithms (sliver-free 
refinement; handling boundary encroachment) to develop our advancing front algorithm. We 
have also extended the analysis by Si m to explain why our algorithm is likely to provide smaller 
high-quality meshes. We will extend this algorithm for constrained Delaunay refinement and 
implement it to study its performance in practice. Other future directions for research include 
adapting the algorithm for surface meshing algorithms and anisotropic meshing in 2D. 
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A Extension to a General Class of Piecewise Linear Complexes 


In our ten-page abstract above, we assumed that the local sizing field at every input vertex is 
equal to (or smaller than) the length of the shortest edge in the input PLC. There are PLCs 
for which the condition does not hold. For instance, consider a line segment ab and a point 
p close to it. If no vertices other than a and b are closer to p, the local sizing field at p is 
equal to its distance to ab, which may be smaller than pa and pb. We will first analyze why 
our algorithm may not function as an advancing-front algorithm in such instances and then 
provide a suitable preprocessing algorithm to achieve our objective of placing vertices in an 
advancing-front manner. 

Let us assume that Aabp is a poor-quality triangle. If we were to carry out Steiner vertex 
placement as in our previous algorithm, Aabp would be split in an appropriate order depending 
on \pa\ and \pb\. Further, if the Steiner vertex for Aabp encroaches upon segment ab, a new 
vertex m would be added on ab, and all free vertices in the circle/sphere with m as the center and 
having a radius |ma| (or \mb\) would be deleted. Since p is not a free vertex, it is not deleted even 
if it is inside the circle. Thus, we would create a triangle with a smaller edge length \mp\. This 
sabotages our idea of introducing progressively longer edges after each Steiner vertex insertion. 

In order to tackle the issue above, we develop an algorithm to preprocess the PLC to add 
auxiliary Steiner points into the PLC such that the diametral/equatorial circles/spheres on the 
boundary edges and facets do not contain the input vertex. Thus, whenever a new vertex is 
introduced, the length of the shortest edge of a new triangle/tetrahedron can only be greater 
than the length of the shortest existing triangle/tetrahedron. Such a preprocessing technique 
would result in an advancing-front Steiner vertex insertion routine. Since additional vertices 
are added into the PLC, the constant associated with the size optimality will be compromised. 

Let us first consider the 2D problem (see Fig. [6]). If an input vertex p lies inside the diametral 
circle on line segment ab, the perpendicular projection of p on the line joining a and b lies on 
the line segment joining a and b. Let the point of projection be m. Further, Zapb < 7r/2. If 
the point m is added to the PLC, p does not lie inside the diametral circle of edge am or mb, 
but m may not be the most optimal location of an auxiliary Steiner point because it may be 
arbitrarily close to a or b. Thus, we perturb the vertex m to a location m such that p does not 
lie inside the new diametral circles, and its location maximizes the minimum distance from all 
other vertices. 

Lemma A.l. In 2D refinement, rn may be placed such that the constant associated with the 

size optimality does not reduee, i.e., \am \ > \pm\ andp does not eneroaeh upon the line segment 

/ 

am . 

Proof. If p encroaches upon the segment am', \am'\ > 2\pm\ and m must lie inside the line 
segment am . Thus, m may be placed such that \pm\ < \am \ > 2|pm|, and we obtain a 
triangulation where p does not encroach on the boundary segment am . If we process vertices 
in the increasing order of their local sizing field, no other input vertex will be closer than the 
distance \pm\ to m . Thus, in 2D refinement, the constant associated with the size optimality 
does not reduce. □ 

We recommend the same algorithm for 3D PLCs as well, but we do not yet have a bound 
on the size optimality in one of the cases described below. As above, we process the vertices in 
the order of their distance to the next nearest input features, i.e., the local sizing field at the 
vertex. Consider a point p that is encroaching upon the triangular facet Aabc. Let the point 
of projection from p onto the plane formed by Aabc be m, and let us assume that m G Aabc. 
The new vertex may be added at m if it is a sufficient distance from other vertices. Since p is 
right above m, no circumsphere of any triangle containing the vertex at m can encroach upon 
it. Further, due to the Delaunay triangulation of the facet, m does not lie not inside any of the 
circumcircles of other triangles. If m is arbitrarily close to one of the vertices of the triangle. 
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(a) input (b) Steiner point (c) input (d) Steiner point 

Figure 6: The placement of auxiliary Steiner points if the local sizing field is smaller than the 
length of the shortest edge in the initial Delaunay triangulation. Two cases are considered here. 
First, vertex m is placed at the point of projection as it sufficiently far from other vertices and 
the diametral circle in the resulting triangulation does not enclose input vertex p. Second, m 
may be arbitrarily close to vertex a, so it is perturbed to m so that we eventually obtain a 
size-optimal mesh. 


say, a, we perturb m to a new location m on the line joining a and m. We will show that |am | 
can be as large as 2\pm\ before p lies inside the equatorial spheres of the triangular facets. If m 
is too close to some other vertices, some other point on the line segment mm may be chosen 
such that it maximizes the minimum distance from all other vertices. 

Lemma A.2. If\am'\ < 2|pm|, no diametral sphere of triangular faeets eneroaehes upon vertex 

p. 

Proof. Consider the unique plane perpendicular to the plane of Aabc and on the line joining 
points a and m. Consider the cross-section of a equatorial (hemi)sphere on the plane. This cross- 
section is similar to Fig. [6] (c) and (d) if diametral semicircles are added to it. By construction, 

I am I > |om|. If p is inside an equatorial sphere, the radius of the sphere is greater than or 
equal to \pm\ because the center of the sphere lies on the plane of Aabc. The cross-section of the 
(hemi)sphere forms an arc smaller than or equal to a semicircle with the radius less than \pm\/2. 
Note that no arc on the plane may enclose vertices a or m because the triangular facets form 
the Delaunay triangulation. Thus, if |om | < 2\pm\, no equatorial sphere of triangular facets 
encroaches upon vertex p. □ 

We still have not considered the following three cases: (a) when point m lies outside all the 
triangles whose equatorial spheres enclose p, (b) vertex p lies close to an input line segment, 
and (c) point m or m lies very close to an input line segment. Case (a) cannot happen 
because the vertex p would have to be inside an equatorial sphere of only one of the two such 
intersecting spheres, but the circle of their intersection projects onto the common edge of the 
two triangular facets. Thus, if the vertex p is inside the equatorial sphere of one of the triangles, 
but its projection is outside, it is also inside the equatorial sphere of the triangle over which its 
projection falls. For case (b), we may follow the algorithm above used for 2D refinement. We 
will provide some nonrigorous analysis for case (c) below. 

A natural solution for case (c) is to place the vertex m at a point on the input line segment 
furthest away from existing vertices such that p does not lie inside the resulting equatorial 
spheres after the insertion of the new vertex. What we do not yet know is how close to an 
existing vertex (relative to the local sizing field) such a point can lie. Consider the unique 
plane that intersects the encroached upon triangular facet on the input line segment and passes 
through the vertex p. Since we know that m or m is closer to the input line segment than the 
distance \pm\, we know that the dihedral angle between our new plane and the plane of the 
triangular facet is less than 45 degrees. Thus, the cross-section of the equatorial (hemi)sphere 
(of a triangular facet incident on the input line segment) on the new plane forms a petal. If 
p is inside the petal, the radius of the circle from which the petal is formed is constrained 
(but smaller than 2\pm\). The constant associated with the constraint results in a proportional 
change in the constant associated with the size optimality of the mesh. 
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We have not yet considered other equatorial spheres that may enclose p after m" is inserted 
into the mesh. If p is still inside an equatorial sphere, it cannot be the equatorial sphere 
associated with a triangle that is incident on the input segment considered above (because 
of our choice of m ). From the reasoning for case (a) above, since p does not lie inside the 
equatorial sphere on the triangular facet in which its its projection lies, it does not lie in the 
equatorial sphere of any other triangular facet on the same plane. The analysis is nonrigorous, 
and we leave the rigorous analysis as future work. 

B Angle Maximization-Based Vertex Placement 

One of the disadvantages of the distance maximization-based routine provided in the paper is 
that the definition of sliver is necessary to place the vertex. Here, we present a slightly different 
algorithm where we attempt to improve the dihedral angles between the facets of a tetrahedra 
by maximizing the weighted distance from the potential new tetrahedra that may be formed as 
a result of the insertion of the new vertex. 

For the angle maximization-based routine, we maximize the minimum weighted distance 
from the plane of those triangles that may form small slivers. The weights are proportional to 
the area of the triangles. Since the volume of the potential sliver is proportional to the area of 
the triangle and the distance of the Steiner vertex from the plane of the triangle, the increase 
in the weighted distance also increases the dihedral angle as well as the volume-edge ratio. As 
mentioned in Section 2, there may be, at most, a finite number of such facets. Li [7] showed 
that a small perturbation of the Steiner vertex may change the connectivity in the resulting 
Delaunay tetrahedralization. Thus, it is important to consider all potential triangular facets 
that may form small slivers. 

For single-vertex insertion, however, this technique can place points that are too close to 
the existing vertices. Hence, the additional geometrical constraint in the form of the picking 
region is necessary. Our algorithm maximizes the weighted distance within the intersection 
of the picking region and the snow globe to insert the Steiner vertex. In the case of vertex 
encroachment, the same modification as above is applied. Although we expect better dihedral 
angles from this routine, it is likely to produce larger meshes. 

For multiple-vertex insertion, instead of the picking region constraint, we place another 
constraint that does not allow vertices to get too close to each other. Specifically, we do not 
allow two vertices to be closer than a/max, where a is some feasible constant and /max is the 
larger of the shortest sides of the two poor-quality tetrahedra that the vertices are attempting 
to eliminate. This constraint is necessary to guarantee the termination and size optimality of 
the algorithm. 

C Additional Analysis 

The ordering of petals and snow globes is the reason why the Delaunay refinement algorithms 
provide smaller meshes. The algorithms attempt to get rid of poor-quality elements with short 
edges as quickly as possible. 

Lemma C.l. Assuming boundary encroachment does not play a role, for petal- and snow globe- 
based refinement algorithms, poor-quality triangles or tetrahedra assoeiated with the shortest 
edge in an existing poor-quality element ean be eliminated with a finite number of Steiner vertex 
insertions. 

Proof. For 2D refinement, with at most two additional Steiner vertices, poor-quality triangles 
associated with a short edge may be eliminated by placing the vertices on the either side of the 
shortest edge within the respective petals and the picking regions. For 3D refinement, because 
we eliminate slivers by avoiding the forbidden regions, at least one new tetrahedron is of good 
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quality with bounded dihedral angles. Thus, after a finite number of Steiner vertex insertions, 
all poor-quality tetrahedra “around” the short edge are eliminated. □ 

If circumcenter insertion is used, there is no guarantee that the new triangle/tetrahedron 
is of acceptable quality, and the triangle/tetrahedron may have to be refined again. We prove 
the following lemma only for 2D refinement. The existence of slivers makes it hard to prove the 
lemma for 3D refinement. 

Lemma C.2. Assuming boundary encroachment does not play a role, for circumcenter-based 
refinement algorithms, poor-quality triangles associated with the shortest edge in an existing 
poor-quality element can be eliminated with a 0(log (D//min)) Steiner vertex insertions, where 
D is the diameter of the domain and Zmin is the length of the edge. 

Proof. Let us assume the radius-edge ratio of the poor-quality element is po- After the first 
vertex insertion, the quality of the new element of the shortest edge becomes pi. After the 
second vertex insertion, it becomes p 2 , and so on until m vertices. We know that pk/pk+i < 2, 
0 < k < m. Since we assume boundary encroachment does not play a role, po is bounded by 
0(T)/Zmin)- Thus, m = 0(log (D//min))- □ 

Theorem C.3. The size of the mesh is 0(log (D/lsfmin))? where D is the diameter of the 
domain and Isfmin is the minimum local sizing field in the domain. 

Proof. Asymtotically, the size of the mesh does not change due to our algorithm. Thus, the 
proof is identical to the one provided by Si m- The constant associated with the size optimality 
is different (as show in the ten-page abstract) due to the ordering. □ 
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